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Abstract 

In this paper we analyze inviscid aerodynamic shape optimization problems gov- 
erned by the lull potential and the Euler equations in two and three dimensions. The 
analysis indicates that, minimization of pressure dependent cost functions results in 
Hessians whose eigenvalue distributions are identical for the full potential and the Eu- 
ler equations. However, the optimization problems in two and three dimensions are 
inherently different. While the two dimensional optimization problems are well-posed, 
the three dimensional ones are ill-posed. Oscillations in the shape up to the smallest 
scale allowed by the design space can develop in the direction perpendicular to the flow, 
implying that a regularization is required. A natural choice of 6ucli a regularization is 
derived. The analysis also gives an estimate of the Hessian's condition number which 
implies that the problems at hand are ill-conditioned. Infinite dimensional approxima- 
tions for the Hessians are constructed and preconditioners for gradient based methods 
are derived from these approximate Hessians. 


‘This research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NASl-19480 while the authors were in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), Mail Stop 13 2C, NASA Langle; Research Center. Hampton. VA 23681- 
0001 . 
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1 Introduction 

In recent years there has been a growing interest in solving optimization problems governed 
by the Euler and the Navier Stokes equations [lM9j. The new interest in this classical field 
[10, 1 1] is due to the increase in computer’s speed and improvements in algorithms for the 
numerical solution of the flow equations. 

The problem of designing a three-dimensional wing requires solving an optimization prob- 
lem with many design parameters. Such a problem may be computationally difficult depend- 
ing on the cost function's level curves in the vicinity of the minimum. A measure for the 
level of difficulty is the condition number of the Hessian. The eigenvalues of the Hessian 
(which is a symmetric operator) are the curvatures of the cost function in the principal di- 
rections. A large deviation in the eigenvalues means that the cost function has level curves 
which are thin ellipses. This is well known in the optimization literature to cause slowness 
of convergence toward the minimum for gradient-based methods [12j. 

Aerodynamic optimization problems are ill-conditioned as noted in [6. 9j, and as will 
be shown in this paper. Therefore gradient descent methods will be extremely inefficient 
especially when the number of design variables is large. A standard method to overcome this 
difficulty is the New’ton method where the Hessian is computed explicitly [12]. The Newton 
search direction is the gradient multiplied by the Hessian's inverse, a computation which is 
impractical in aerodynamic optimization problems since it involves numerous solutions of 
the flow PDEs. On the other hand jsing low rank quasi-Newton methods, such as BFGS 
[3, 4). will result in a deteriorate con -ergence as the number of design variables increases. 
Thus, a new method is required. 

Another difficulty in inviscid aerodynan,’ optimization problems is the ill-posedness 
of three dimensional problems which shows up as small scale oscillations in the shape in 
the direction perpendicular to the flow. Such oscillations were observed in applications as 
reported in (9). One way to avoid these oscillations is to apply smooth finite dimensional 
representation of the shape in the spanwise direction. Another approach is regularization 
by introduction of a penalty to the cost function for oscillations in that direction. The 
need for penalising the cost function in order to remove oscillations was observed by [3] 
However, iri that case the oscillations were a result of the discretization and had no differential 
counterpart. Penalization was used also in two dimensions where the differential optimization 
problem is well-posed. 

In this paper we develop a new approach to approximate the Hessian and its inverse for 
optimization problems governed by PDEs. Hessian symbols were previously computed fer 
smoothing predictions in the development of multigrid one-shot methods [13]-[16]. Here, a 
similar analysis is applied to inviscid flow problems including the full potential and Euler 
equations in two and in three space dimensions. In Sec. 2 the optimization problem is defined 
together with its small disturbance approximation. The necessary conditions for a minimum 
and their relation with the Hessian are discussed also in the finite dimensional design space. 
In Sec. 3 local mode analysis is presented to approximate the Hessian's symbol. The analysis 
is local and. involves freezing the coefficients to obtain a problem in half space with constant 
coefficients, where Fourier techniques are employed. In Sec. 4 the analysis is applied to an 
optimal shape design problem governed by the full potential equation. In Sec. 5 the analysis 
is applied to the Euler equations, and a symbol is obtained identical to the full potential 
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case. In Sec. 6 the Hessian's symbol is analyzed and conclusions are made concerning the ill- 
conditioning and ill-posedness of the problems at hand. In Sec. 7 regularization is discussed 
to avoid the ill-posedness which exists in three dimensions in the spanwise direction. In Sec.8 
preconditioners are developed for the small disturbance and the optimal shape problems in 
subsonic and supersonic flow. Finally in Sec. 9 discussion and concluding remarks are made. 

2 The Optimization Problem 

A typical inviscid aerodynamic optimal shape design problem aims ai finding the shape of 
a surface, e.g.. airfoil or wing, such that the resulting pressure distribution on that surface 
trill minimize the least squares distance from a prescribed pressure distribution. Let 17 be 
a domain in Ft? and T(;r) a parametric representation of the part of the boundary dQ to 
be designed. The optimal shape problem is to compute the boundary position, T. that 
minimizes a cost function defined on F, e.g., 

min J^f(U)do (2.J.) 

where / is a prescribed function and U is the solution of a PDE defined on ft, 

C(T,U) = 0. 


2.1 The Small Disturbance Approximation 

For the analysis of the Hessian it is enough to consider small perturbations of the boundary 
T. In order to further simplify the derivation we consider a localization of the problem in a 
vicinity of a boundary point and study the resulting half space problem. Let us introduce 
the following notation 

E” + ={(.r,x n ):r<zlR n - t . 0 < * H } 
s {(x..r„) :f € ZFT” ] ,x n = 0}. 

We consider perturbations of the form 

r - F* + ean 
U = U' + eO -f 0(f 2 ) 

where F“ and U" are the optimal boundary shape and state solution respectively; n is the 
outward normal and e is a sufficiently small positive number. The resulting optimal control 
problem is obtained using a Taylor expansion in c and is discussed in more details in the 
next sections. The new control problem, for the control variable at, is defined in half space 
and is referred as the “small disturbance minimization problem’*, namely, 

min F(q, f ; ) (2.5) 

subject to the following state equation : 

L V ~ 0 on JR* (2.6) 

B V = Cd on (2.7) 


(2.2) 

(2.3) 

(2.4) 
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and additional homogeneous conditions on the solution at infinity. For this control problem 
the minimum is attained for d = 0. 

The necessary conditions for a minimum of the small disturbance problem are given 
by the state equation Eq.(2.6-2.7) and two additional equations: design and costate. The 
design equation is a PDE defined on the boundary dlR\ and is denoted by, 

.4{£\A,a) = 0 on dR\, (2.S) 

where V is the solution of the state equation (2.6-2. 7) and A is the solution of the costate- 
equation . i.e., 

LA = 0 on 1R Z + 

BA = tO on dM\ 

and additional homogeneous conditions at infinity. 

We assume that the perturbation d is composed of high frequencies and compute the 
resulting solutions, V and A, in the vicinity of some arbitrary point on the boundary <9iR|. 
The solution there is approximated by a constant coefficient, problem defined in half space 
where Fourier analysi can be applied. 

It can be shown that for feasible solutions of the state and costate equations. U = t'’(d) 
and A — A(t~’(d)).. the design equation residuals are equal to the gradient of the cost; function 
with respect to the design variables: 

g(6j'Ha)) = -4(t/(d), A(£(d)),d). (2.11) 

In the vicinity Tie minimum the following relation holds 

g(a.C'(a)) - ^(0,0) + Ha -f h.o.t. (2. 12) 

where we denote by H the Hessian, i.e., 

II = Vzg{ 0,0). 

A Taylor expansion, of the right hand side in (2.11) and comparison with (2.12) results in 

H = AqO(, + A\A{jO& + A*. (2.13) 

The dimension of the ilessian. H. is determined by the dimension of the design variable, a. 
If d belong to a finite dimensional space of dimension N then the gradient is a vector of size 
N and the Hessian is a N x A 7 matrix. If d belongs to an infinite dimensional space (e.g., 
some function space) then the gradient is an element in an infinite dimensional space and 
the Hessian is an infinite dimensional operator. The Newton step, a, satisfies 

H6 = -g (2.14) 

where g is the gradient at the given iteration. 


(2-9) 

( 2 . 10 ) 
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2.2 Finite Dimensional Design Space 

In applications it is a common practice to restrict the design space to a finite dimensional 
subspace where the shape F. or in the small disturbance model o, is given as 

(2.15) 

where N is the dimension of the design space, c*j are the design variables and f } are fixed 
basis functions. The gradient for the finite dimensional case, g — (<?j, • • • , g.v), is obtained 
by a projection of the infinite dimensional one onto on the finite dimensional space, i.e., 

Jg = b (2.16) 


where J is a matrix defined by 


Jii = (/.(»). U = (2.17) 

and the right hand side, 6, is a vector whose elements are given by 

M 8 ) = (»(«), /*( 8 ))i* k - 1. • • • , A'. 

Let hjk be an element in the matrix representing the Hessian using the basis functions 

hjk =Wj. /*)*»■ (2.18) 

A Newton step for minimizing the cost function consists of moving in the direction a = 
(ttj, •• *,a,v) given by: 

ha = -b. (2.19) 

Note that this equation is obtained by projecting the general equation (2.34) onto a finite 
dimensional space spanned by {/jJjlj. 


3 Local Mode Analysis 


A local analysis of the operator H is done next. The Fourier analysis consists of analyzing 
the solution of the following system of equations in half space: 


LU — 0 on 81+ (3.1) 

B 0 = C a on dRl (3.2) 

La = 0 on Rl (3.3) 

BA = £0 on dR\ (3.4) 

A(U,&.a) = 0 ondJR n + , (3.5) 
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where we replaced 0 by 0 to account for cases in which L stands for a system of PDEs. 

The computation of the symbol of the Hessian. (2.13), near the minimum is done by 
considering a perturbation in the design variable of the form 

a(f) = 6 (k)e'** (3.6) 

and calculating the corresponding term Ha. The small disturbance solution U can be rep- 
resented as 

h ?. *») = £ p J (k)i / j (k)e rkf e ikix \ (3.7) 

where q equals the mtniber of boundary conditions in (3.2). Each of the terms in the 
expression for U satisfies the equation 

LV ? ,-(ib)e it, e tti ** « 0. (3.8) 

which implies 

L(k.ki)i' ] (k) = 0. (3-9) 

where t(k. &£) is the symbol of L. Moreover, Eq.(3.9) implies that for j ~ !•••<?, 

6etL{k,k{) - 0 (3.10) 

and that Vj(k) is an eigenvector of L(k, k{) witli a zero eigenvalue. Substituting the expres- 
sion for V into the boundary condition (3.2) results in the following equation. 

£>,(*)£(*, ^)?,-(*) = C{k)a(h (3-11) 

K -A 

where B(k. k J „) is the symbol of B. and C(k) is the symbol of C. Introducing the matrix 

W(k) = (B(k. kl)t\(k), .... B(k, *')?,(£» 


and the vector 



Eq.(3.11) can be written as 




W(k)S(k) = C(k)a(k), 

(3.12) 

or, equivalently, as 

U(k) = V(k)W~'(k)C(k)*{k). 

(3.13) 

where 0(k) is defined by 

0(2. 0) « V(k)8(k)e ilx e 0(k)e ili 
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and where V‘(£) is the matrix V' = (V\ (fc), . . . , V%(£)). 


The adjoint equations are treated in a similar way. The solution A is represented by 

A(£, *,) = £ 3j(k)kj(k)e^' r t' i ' ,Ir ' (3.14) 

i=t 

where q + q — N and N is the degree of the polynomial in (3.10). Let k } n be the roots of 

detL(k,k{) = 0 (3.15) 

and z.j the eigenvectors of L(k. ££) corresponding to a zero eigenvalue. The costate boundary 
condition (3.4) implies 




(3.16) 


•=i 


Introducing the matrix 

mb = <M. £)!,(*), . . . , ^(fr. **)!,(*)) 

and the vector 

m = ( A( ?)= — js # (*)) 

Eq.(3.16) can be written as 

#(*)l(*) = £(&)£(*), (3.17) 

or. equivalently, as 

•*(*) = t(k)W~ 1 (k)C(k)C 7 (k). (3.18) 

where E(£) is the matrix ~.(k) = (Ej(k), Ef(k)). 


Substitution of (3.13) and (3.18) in the symbol of the Hessian’s expression (2.13) results 
in the following formula for the Hessian’s symbol: 


H(k) = A 0 YW- l C + Aa 


w^Cyw-'c + Az. 


(3.19) 


4 The Full Potential Equation 

In this section we apply the ideas discussed above to minimization problems governed by the 
full potential equation. We consider the problem of optimizing the shape of an aerodynamics 
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configuration so that the shape minimizes the deviation of the pressure from a target pressure 
distribution, namely. 

min F(T,p) - j^(p-p“) 7 dcr (4.1) 

subject to 


£(P. o) = div(pV$) = 0 on 0 
d>„ = 0 on F 


(4.2) 

(4.3) 


and an additional boundary condition at the far-field (the notation stands for the outer 
normal derivative of 0 on the boundary). 

The density, p. is given by the isentropic density law [17;: 


ajj. 

poo * 2 H x J 


(4-4) 


with 7 the specific heat ratio and p^ and H^, the values of the density and total enthalpy 
at infinity. 

The pressure, p, is related to the density p and the speed of sound, c, by 





(4.5) 


where for perfect gases the speed of sound is related to the potential through the relation 


£»«(?- l) r ff, c~ 


(v<>) ! 


2 1 ' 


(4.6) 


The Mach number M is given by 


M 2 = 




!2 


4.1 The Small Disturbance Minimization Problem 

The derivation of the small disturbance minimization problem follows the argument and 
notation of Sec. 2 . 1 . We perform localization and set the local coordinate system on the 
boundary such that the flow is in the x-direction (Vd> — (<^ r . 0 , 0 )). On the perturbed 
boundary a Taylor expansion gives [ 1 ] 

= ? + p. + ap n ±h.o.i. 

and 

drf ptrtM „ (! _ t)da + k.o.t. 

H 

where = jj- -f and .ftp J ?2 are the principal curvatures. 

The small disturbance minimization problem is giver: by 

mn\F(a,U)~ f ({p - p~) + p(4>) + &p n ) 2 dxdy - f ~{p - p*) 2 dzdy. (4.7) 
» Jd r 3 , Jnn\ R 
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subject to 


(4.8) 


(1 - A/ 2 )<p x - - = 0 on JRi. 

Tbe small disturbance boundary condition is given by 

^ (p^otg ^ (p*xOt s 0 (4*9) 

where the unperturbed flow is in the x-direction. Since we are analyzing the Hessian for the 
high frequencies it is enough to consider 

- 6, = 6 x a x . (4.10) 


The relation between p and <p in the cost function (4.7) is obtained from the relations (4.4- 
4.6). We obtain that the relevant cost function in terms of d to be minimized is 

/( q, <?) ~ I t (( p~p’)~ P4>x<j>x) 2 dxdy - f a [~(p - p" f - 2 p n ( p -p‘)j dxdy. (4.11) 

4.2 The Adjoint Equations 

By standard variational calculus it can be shown that the gradient of the 
given by 

(p ~~ p" ) 2 

g = c> r A r + 2 p n (p - p*) 

where A is the solution of 

(i — A/ 2 )A IX -f- A yy + A. r = 0 on IFtt (4.’!!} 

with the wall boundary condition 

- \ t - 2 [p<M(P ~ P*) ~ P<Px^r)] x = 0 on dlFL\ (4.14) 

and <j> is the solution of (4.8) with the boundary condition (4.10). We also require that the 
solution is bounded for the subsonic case (M < 1) and that in the supersonic case (M > 1) 
no waves propagate in the direction of —Vd{x 0 ), i.e., not propagating in the negative x- 
direction. This requirement is done so that the solution of the small disturbance problem 
will be consistent with the far-field boundary conditions of the unperturbed problem. 


cost function is 
(4.12) 


4.3 Local Mode Analysis 

We now go through the analysis in Sec.3 in order to compute the symbol H{k) (see Eq.3.19) 
using the full potential state equations. Following a perturbation a 

a{x) = d(k)e' ls 

the small disturbance solution, d. can be represented as 

m = MW lt e ,k * s , (4.15) 
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and. similarly, the adjoint variable, 

A(.f) = \(k)e i!: *c ik> \ (4.16) 

Substitution of (4.15) in (4.8) and of (4.16) in (4.13) results in the following relations for 
k-> and for i 3 : 

( 1 - M 7 )k'l + kl -f A'| = 0 (4.17) 

( 1 - M 2 )** + k l +k 2 = 0. (4.18) 

These are the analogs of equations (3.10) and (3.15) respectively. 


Y 


incorrmp 


t 

! 



X 


outfoing 


Figure 1: The supersonic flow can be decomposed into two waves: incoming and outgoing 
the plane 

The choice of sign for k 3 and k 3 is done as follows. Since the half space problem is related 
to localization of the original problem around some point xo € T, the solutions that we 
construct in half space should be compatible with the far field boundary conditions of the 
original problem. In the supersonic regime the solution can be decomposed into two waves: 
incoming and outgoing of the plane (see Fig.l). In terms of the local coordinate system, 
the incoming characteristic has a component, in the negative y-direction while the outgoing 
has a component in the positive y-direction (both propagate in the positive x-direction). 
Since perturbations of the shape F can not change the far field inflow data, the change in the 
incoming characteristic should vanish. Thus, the part of 4> which propagates to the negative 
x-direction is set to zero. This implies that k 3 is of opposite sign to that of k\ when k 3 is 
real valued (i.e. supersonic flow). The. adjoint variable has characteristics in the opposite 
direction and therefore we require that the part of the solution for A which propagates to 
the positive x-direction should be set to zero. As a result the sign of k 3 is the opposite of 
that of k 3 if they are real valued. In the subsonic regime k 3 and k 3 are imaginary therefore 
the proper sign is positive for both so that exponentially decaying solutions are obtained. 
Therefore k 3 k 3 = — 1& 3 | 2 for both subsonic and supersonic flow. 

From the boundary condition (4.10) we obtain a relation between a{k) and <£(fc), (as in 
(3.13)) 


m = 

"3 
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(4.20) 


and from the boundary condition (4.14) (as in (3.18)) 


\(k) = 2ip 2 <i>l&6(b 

"3 

Hence, the gradient (4.12) in Fourier space is 

g(k) = <t> x iki\(k) -i- l.o.i. 

Substitution of (4.19) in (4.20) and of (4.20) in (4.21) results in the relation 

g(k) - 2 p 7 0 A x ~ti - ~2p 2 <j> 4 kl ~ 




O' 


(4.21) 


(-4.22) 


hh ]At-3 

and thus by relation (4.17) we obtain the following formula for the symbol of the Hessian: 


H(k) = 2 p 2 o A x 


|(1 - M*)k\ + fcff 


(4.23) 


5 The Euler Equations 


As in the full potential case we consider the problem of optimizing the shape of an aerody- 
namic configuration, subject to the Euler equations, such t hat the optimal shape minimizes 
the deviation of the pressure from a target pressure distribution (4.1). We perform the anal- 
ysis away from shocks so that it can be done using a non-conservative formulation. The 
Euler equations in quasi-linear non-conservative form are given by 


C(T,U) = 


Q pd x pd y pd, 

0 Q 0 0 

0 0 Q 0 

0 0 0 Q 

V 0 pc 2 d x pc 2 d y pc 2 d t 


0 

A 

l A 

l A 

Q 


\ 



u 


V 


1 0 

/ 

\ p J 


o 


(5.1) 


where Q = u ■ V (» ~ (u. v,iv) denotes the velocity vector), with the solid wall boundary 
condition 


u-n = 0. (5.2) 

Additional conditions that are given at the inflow and outflow boundaries in terms of char- 
acteristic variables are not used explicitly in the derivation of the approximate Hessian. 

5.1 The Small Disturbance Minimization Problem 

Following the same argument given in Sec.4, the small disturbance cost function is given 
in Eq.(4.7). The stale equations (5.1) are perturbed in the vicinity of the minimum and 
t<; » perturbation variables solve the linearized Euler equations which up to low order terms 
are given by the same matrix operator as above. Following the localization and half space 
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approximation we set the local coordinate system on the boundary such that the flow is in 
the x-direction (u = (h. 0.0)). The small disturbance wall boundary condition is given by 
(the perturbation variables are denoted by p„ v, and f>) 

- tr = udj (o.3) 

where we have omitted as before the zero order terms in ti. 

5.2 The Adjoint Equations 

By standard variational calculus the gradient of the cost function is given by 


g = — u/>(A] -f c 2 A 5 ) x - + 2p„(p - p’) (5.4) 

where ( A] . A. A$) is the solution of the following system of equations in 1R\. 

djr(Aiu)=0 (5.5) 

grad{p\\) — grad{u ■ A) -fc^rcdipAj) = 0 (5.6) 

dir(-A) + div(\$u) = 0. (5.7) 

P 

The boundary condition on dffti is given by 

- -A 4 + 2 p(p -/>') = 0. (5.S) 

P 


An additional requirement is consistency with the far-field boundary conditions of the un- 
perturbed problem; i.e.. the solution is bounded for the subsonic case (M < 1) and that no 
waves propagate in the direction of — u(io) in the supersonic case (M >1). 

Neglecting zero order terms in Eqs.(5.5-5.7) w r e arrive at the following quasi-linear form: 


Aj A 


Q 0 0 0 0 A 


A, ^ 

A 2 ] 


pd x Q 0 0 pc 2 d r 


x 2 

A3 

“ — 

pd v 0 Q 0 pc 7 d v 


^3 

Ai 


pd j 0 0 Q pc*d z 



a 5 ) 


0 A l A l A Q J 


t As j 


5.3 Local Mode Analysis 

Following the same procedure as in the full potential case we consider a perturbation in the 
design variable of the form 

d(x) = «(k)c i£ * (5.9) 

with k = and .r = (x,p). As a result the states and costates are perturbed by 

C(x) = £'(fc)e ,X V^ (5.10) 

A(x) = K(k)^ r e iht (5.11) 
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with U — ip.u.p) and A = (A!, A, A 5 ). Introducing the notation 

<?(£) = u - k = uk x 

and 

/ Q pk i pit 2 pk 3 0 ^ 

0 Q 0 0 ^ 

L{k.h) = i 0 0 Q 0 ^ 

0 0 0 Q & 

£ 

\ 0 pc 2 k\ pc 2 k 2 p<?k% Q / 

the form of the above solutions for the perturbation variables imply that 

deti(A\ k 3 ) — 0 (5-12) 

dett(k,k 3 ) = 0 (5.13) 

where L is the adjoint of L. The two relations result in a fifth-order polynomial equation for 
the wave numbers (k x . k 2 . k 3 ). 

k x 3 u 3 {k\ 2 u 2 - c 2 (/v 2 + k\ -f kl fj = 0 
k^u^ik^u 2 — c l {k\ 4- k 2 + ~ 0, 

with the roots 

k\ — k\ — — (1 — M 2 )k 2 — k\ and k x = 0, (3,14) 

where the Mach number is given by M 2 — Note that the roots for k 3 and for Z- 3 are 
identical to those obtained in the full potential case. 

Let 

V\(k) - (k 1 puc~ 2 .-k u -k 2 ,-k 3 .k l pu) 

V, ( k) — ( hi puc ~ 2 , — kt . — k 2 , A* 3 , k x pu ) 

f 3 (*) = d,o.o, o,o) < 5 - 15 ) 

\\{k) = (0.1, 0,0,0) 

V 9 (k) = (0.0.-A' 3 .L- 2 .0) 

E ](k) = (O, —<?pk\,— (?k 2 p. — c 2 pk 3 .kiV^j 

r. 2 (k) = (o. -c?pk\, -e 2 k 2 p. c 2 pk 3 , k\ii) 

i 3 (j t) = (-e 2 , 0.0. 0,1) (- 5 - 16 ) 

!(*) = (0,1, 0.0,0) 

t s (k) = (0.0, —k 3 , ^2,0). 
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The vector Vi{k) is the eigenvector of L(k. kz[ h)) |Z(fc, &3(£))j with a zero eigen- 

value, the vector (?*(£) [z 2 (i-) 1 is the eigenvector of L(k.-k.z(k)) |Z(A\ -Jt 3 (^))j with a zero 

eigenvalue aod the vectors V 3 . 4i5 (^) jE3. 4 .$(£)] are the eigenvectors of Z(0, A'3) [Z(0. fc 2 , £3)] 

with a zero eigenvalue. 

The eigenvectors which correspond to k t — 0, (V*3, 4>5 (fc) and Ez^^k)), represent waves 
coming from -oc to 00, in the half space coordinates, parallel to the surface dIR \ in the 
analysis problem, and from 00 to ~oc in the adjoint problem and therefore do not play a 
role in the analysis. 

The eigenvectors which correspond to £3 — kz(k) and £3 = k$(k) are relating changes in 

A. 

the designed surface with the flow field. However, only V j (k) and ~i(k) are consistent with 
the far-held boundary conditions as discussed in Sec. 4.3. 

We look for amplitudes of the vector solutions i'(k) in (5.10) and A(&) in (5.11) consistent 
with the boundary conditions. From the boundary condition of the state equation (5.3) we 
get (see (3.12) together with Vf in (5.15)) 

k 3 ,S(k) — uikia(k). (5.17) 


The boundary condition of the costate equation (5.8) (see (3.17) together with the definition 
of 5j in (5.16)) implies 

c 7 kzp(k) — -2kipu8(k) (5.18) 

and the design equation (5.4) implies 

g(k) = puHktfm -f i.o.i. (5.19) 


Finally, substituting of (5.17) in (5.18) and of (5.18) in (5.19) yields 


(5.20) 


6 Analysis of the Hessian 


In the previous sections the symbol for an approximate Hessian near the minimum was 
obtained, namely, 


H(k\, ki) = 2 p 2 u 4 h(k x ,k 2 ) 


( 6 . 1 ) 


with 


](4 

h(h,h)= + k r { - 


( 6 . 2 ) 


The fact that the same Hessian is obtained both for the full potential equation and for Euler 
equations implies that for the purpose of developing new optimization algorithms it is enough 
to consider the full potential equations. Although the Euler equations presents additional 
difficulties compared with the full potential equation for the analysis problem no additional 
ones exist as far as the optimization is concerned. 


13 


ADA309055 



The Hessian’s symbol in the discrete space can be obtained in an analog way. The result 
will then depend on the specific discretization we use to solve the equations. One can get 
a crude approximation of the discrete Hessian by replacing the wave number fcj with the 
discrete wave number Oi (and k 2 with 9 2 ). 


(*„* 2 ) = (£-, 

hi 



where {h\.h 2 ) are the mesh-sizes in the (x.y) directions respectively. 


6.1 Two Dimensions 


In two dimensions the x 2 direction does not exist and thus in the Fourier space we set k 2 = 0 
in Eq.(6.2) resulting in 


h(k j) = 


*? 


|1- iVP|' 


(6.3) 


Substitution of k j = in (6.3) implies that the condition number of the Hessian 3eales with 
the grid mesh-size as 0(-^). Thus, the Hessian is ill-conditioned and its condition number 
increases quadraticallv in the discretization parameter as the grid is refined. Therefore, a 
good estimate of the Hessian is required to obtain a fast convergence of the optimization 
process. 


6.2 Three Dimensions 

In three dimensions the properties of the Hessian are det ermined by 

k 4 

« j^ (1 _ a ’ /2) + JTy 

For fixed wave number in the stream direction, = const, the symbol approaches zero 
as the wave number k 2 approaches infinity. This means that the cost function is nearly 
flat with respect to perturbations in the shape which are highly oscillatory in the direction 
perpendicular to the flow. This might explain recent numerical results showing that the 
wing surface is likely to develop oscillations in the spanwise direction [9] . Note that the 
above oscillations do not appear in the 2D problem which indicates that the 3D aerodynamic 
optimal shape design is inherently a more difficult problem than the 2D problem. 


7 Regularization 

As discussed in Sec.6.2 the three dimensional problem is ill-posed and as a result oscillations 
are expected to appear in the direction perpendicular to that of the flow. In order to 
eliminate this phenomenon it is necessary either to penalize the cost function or alternatively 
to represent the shape as a finite sum 'X base functions which are smooth in the direction 
perpendicular to that of the flow (roughly the spanwise direction). 
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In this section we propose a natural penalty of the cost function which will prevent the 
oscillat ions withou' increasing the computational cost of the preconditioner. The penalty 
term has the following form. 

F r ~ 7 j r (Tc*) 2 da (7.1) 

where rj is a positive parameter and T is an operator whose symbol is 


fikukt) 


k\ 


This results in a modified Hessian given by 


H{k l ,h)±2r J 


g 

\kt(\ -M*) + P 2 \ 


(7.2) 


which reaches asymptotically a constant as k 2 increases, for fixed value of k\ . Thus, for large 
values of |&i | 4 - |A* 2 | we have > 6 > 0 for some positive 6 , which implies that 

for the modified problem the high frequency perturbations in the shape are well behaved. 
As the shape f converges to the minimum, 7 can be decreased, resulting in a solution 
which is regularized in the direction perpendicular to the flow. The implementation of the 
regularization is discussed in the next section. 


8 Preconditioning 

In Sec.6 we concluded that aerodynamic optimization are ill-conditioned problems. Therefore 
having second order information is crucial for effective convergence. However, the explicit 
solution of a Newton step. H& — —g, requires to compute the Hessian, H, explicitly and then 
to invert it. This will become practically impossible for a realistic aerodynamic optimization 
problem computed numerically on a fine mesh and having a. large number of design variables. 
Using low rank quasi-Newton methods, such as BFGS. will deteriorate as the number of 
design variables increases. We therefore suggest to approximate the Newton step in the 
differential level (infinite dimension), using the Hessian’s symbol, and then to project the 
result onto the finite dimensional design space which is used in practice. 

The equation defining the Newton direction, for a. in Fourier space is given by 

or explicitly by using the symbol (6.1) 

(2pV*< 4 27*f)d(* 1 , k 2 ) = -\kftl - M 2 ) 4 k 2 \g(k„ k 2 ). ($.1) 

The symbol in the right hand side of equation (8.1) corresponds to a non-local operator in 
the real space. The term which multiplies 7 accounts for the regularization penalty term 
discussed in the previous section (7 should be set to zero in two dimensional problems' 


16 


ADA309055 



8.1 Preconditioner for the Small Disturbance Problem 

Using the relation bet ween fcj, k x and k 2 , 

|*?(1 -Af*) + 'ji = ~{ih)iih), 


we arriv e at 

(2 A 4 *! 1 + 27)k 2 7 }a(ki, fc 2 ) = [ikiWhWki.ki). (8.2) 


This implies the following equation in real space. 


2p 2 v 


4 &6 

dr* 


n 

\ w; - 1 " 


Hz 


+ M 


on 


dltf+. 


(8.3) 


The terms which multiply p are added both to ensure a unique solution to Eq.(8.3) and 
to account for the low-frequencies. Note that Eq.(8.2) is a good approximation for the 
symbol of the Newton equation only in the high-frequencies. In the low- frequencies the 
terms multiplying p, in Eq.($.3), are dominant and result in a steepest descent step, while in 
the high-frequency regime they diminish and a Quasi-Newton step is taken. The term t 
satisfies the following coupled PDE system 


(I- 

- - 

■ !f + »<;> = o 

on IR^ 

(8-4) 



ii 

^ 5 

on dfit\ 

(8.5) 

(1- 

- - v’ff * <l'iV ~ 0 

on ]R\ 

(8.6) 



*K> 

II 

to 

on . 

(8.7) 


We also require that the solution be bounded and that in the supersonic regime (M > 1) no 
waves propagate in the positive x-direction in Eq.(8.5) and that no waves propagate in the 
negative x-direction in Eq.(8.7). Note that the operator T in Eq.(7.1 ) need not be explicitly 
evaluated. By adding the t] term in (8.2) and solving Eqs.(8.3-8.7) we account for such an 
operator. A similar preconditioner can be derived for the small disturbance Eiuler equations. 


8.1.1 Purely Subsonic Flow 

In a purely subsonic flow it is unnecessary to go through the above procedure since 


|A?( 1 - M 2 ) -e k \ | = kl(l ~ M 2 ) ~ k\. 
Therefore the following PDE should be solved on F only: 




0 2 a 


d 2 g . d 2 g 


^9? ~ ‘^k: ^ - a - m 2 )^ + + m - 




d 2 x & 2 y 


( 8 . 8 ) 
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$.2 Practical Implementation 

in practice equations (8.4- 8.7) are replaced by full potential equations defined on the domain 
ft as shown in the following. Let <p be the solution of t he full potential equation (4.2). We 
define 

= rtfO) 

^♦(2) — 4- f! /‘(2) 

where £<1. We claim that the preconditioning equation for d is given by 

2p 7 (y<p ■ V^d V) 2 d -f /id = V9 on ^ (8 9) 

where e is a unit vector perpendicular to V^> and to n (roughly in the spanwise direction): 
x e = it. The function is the solution of the following coupled PDE system: 

S7p(v (1> )V^ = 0 on ft 

0 (1) = — " r) + <#? cn T (8.10) 

V* 1 ' = on dfl — r 

W (2 >)V<A< 2 > = 0 on ft 

— eg ± -•? on F (8.11) 

t!-^ — on dCl 0 f 

where SQ C / stands for the far field outflow boundary. The solution, d, of Eq.(8.9-8.11) is the 
preconditioned search direction to be used in optimization algorithms. This search direction 
will avoid oscillations in the shape in the direction perpendicular to the flow and will require 
many fewer optimization steps to solve the problem. For the Euler equations an analog of 
Eq$.(8.10-8.11) can be derived using the corresponding small disturbance preconditioner. 

8.3 Implementation in a Finite Dimensional Design Space 

In a finite dimensional subspace we replace q in Eq.(8.9) by (see (2.15)) 

d(£) «£*,.&(*). 
j=i 

We then take the inner product of Eq,(3.9) with jit for == 1 • • • N, resulting in a linear set 
of equations for o j : 

(h + ///)<*= -b + gy, (8.12) 

where I is the unit matrix and fi is a positive parameter, as in Eq.(8.3). h is an N x N 
matrix in which an element h l} is given by 

*«■ « (V(SV • V)V, - 2,(c ■ 
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a is the finite dimensional design variable a — (dj. and an element in the right 

hand side vector, b, is given by 


6 ,= 


l 

£ 



The solution of Eq.(8. 12) is the preconditioned search direction replacing 





9 Conclusions 

In this work new local mode analysis for optimal shape design problems which are governed 
by PDEs was developed. The analysis was applied to aerodynamic shape optimization 
problems governed by the full potential and the Euler equations. The analysis was done in the 
infinite dimensional design space where arbitrary changes in the wing's shape, in the normal 
direction, are allowed during the course of optimization. In this case the Hessian is an infinite 
dimensional operator defined on a space of functions and its eigenvalue distribution served to 
study the well-posedness of the optimization problem as well as for deriving preconditioners 
to accelerate the numerical convergence of gradient based methods. In practice, however, a 
finite dimensional design space is commonly used for which the Hessian is a finite dimensional 
matrix. The application of the infinite dimensional analysis to finite dimension was obtained 
by a simple projection. 

The analysis is local and uses freezing the coefficient to obtain a problem in half space with 
constant coefficients, where Fourier techniques are employed. The eigenvalue distribution 
of the Hessian is analyzed by computing its Fourier symbol. It was shown that for two 
dimensional flow the Hessian is a second order differential operator ^efined on the designed 
boundary. In three dimension the Hessian is a pseudodifferential operator (non-local) and its 
properties are much more complex. For both the full potential and Euler flow the symbols 
of the Hessian are identical. Therefore, the complexity of the optimization problems is the 
same for both, although the analysis problem for the Euler equations is more difficult. 

The symbol of the Hessian implies that the three dimensional problems are ill-posed, and 
arbitrary oscillations in the shape can develop in the direction perpendicular to that of the 
flow (roughly the spamvise direction). This explains recent numerical results showing that 
the wing surface is highly oscillatory in the spanwise direction (9) . A regularization that 
involves smoothing only in that direction was introduced and analyzed. Also note that if 
the problem were to minimize drag, rather then matching the pressure distribution, then 
oscillations in the spanwise direction are not likely to appear since they will increase the 
surface area of the wing resulting in an increase in the drag. 

The explicit form of the symbol of the Hessian also implies that these minimization 
problems are ill-conditioned and their condition number increases quadratically with the 
dimension of the design space. Therefore gradient descent method will be inefficient and 
second-order information, by approximating the Hessian (or its inverse), is essential for fast 
convergence. However, low rank quasi-Newton methods, such as BFGS, will deteriorate 
as the number of design variables increases. New preconditioners which approximate the 
inverse of the Hessiari are proposed. Their numerical implementation will be presented 
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elsewhere. These precondilioners are of low computational complexity for two-dimensional 
flow and for purely subsonic flows in three dimensions. In non-subsonic three-dimensional 
flow the preconditioning involves the solution of the full potential equation twice per each 
optimization step, though a substantial decrease in the number of optimization steps required 
to reach the minimum is anticipated. The preconditioning of the suggested regularization 
is straightforward and requires negligible additional computational work. In case a finite 
dimensional design space is used the preconditioning requires the solution of a linear set of 
equations. 

Applications of similar analysis to aeroelastic optimization are discussed in [18]. 
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